function [dyi,xi, zi, zsa, zlagsa, pos_xi, pos_ylags, pos_stars]=create_regressors_ci(Y,X, lx, ldy, lcsx, lcsdy, indf, lcsz,tau_range, CS_demeaned_ave, countryi, debt, lz, dflz)


[N,Tmax,kx]=size(X);         % dimensions
k=kx;
ktau=size(tau_range,1);

[a,b,kz]=size(indf(tau_range(1),debt)); % just to get dimensions of the threshold variables


% find Tmin and Tmax (for each CS unit)

Ty=[ones(N,1),ones(N,1)*Tmax];
Tx=repmat([ones(N,1),ones(N,1)*Tmax], [1,1,kx]);
Tz=repmat([ones(N,1),ones(N,1)*Tmax], [1,1,kz,ktau]);

Z=NaN(N,Tmax,kz,ktau);
for ii=1:ktau
    tau=tau_range(ii);
    %debt=X(:,:,1);
    pom=indf(tau,debt);
    Z(:,:,:,ii)=pom;
end

t=1:Tmax;

for i=1:N
    % dependent variable
    Ty(i,1)=max(Ty(i,1), min( t(not(isnan(Y(i,:))))) );
    Ty(i,2)=min(Ty(i,2), max( t(not(isnan(Y(i,:))))) );
    % regressors
    for s=1:kx
            Tx(i,1,s)=max(Tx(i,1,s), min( t(not(isnan(X(i,:,s))))) );
            Tx(i,2,s)=min(Tx(i,2,s), max( t(not(isnan(X(i,:,s))))) );
    end %s
    for sz=1:kz
        for s=1:ktau
            Tz(i,1,sz,s)=max(Tz(i,1,sz,s), min( t(not(isnan(Z(i,:,sz,s))))) );
            Tz(i,2,sz,s)=min(Tz(i,2,sz,s), max( t(not(isnan(Z(i,:,sz,s))))) );
        end %s
    end %sz
end %i

% de-mean data for computation of star variables
Yd=NaN(size(Y));
Xd=NaN(size(X));
Zd=NaN(size(Z));
if CS_demeaned_ave==1
    for i=1:N
        yp=Y(i,Ty(i,1):Ty(i,2));
        tp=Ty(i,2)-Ty(i,1)+1;
        Yd(i,Ty(i,1):Ty(i,2))=yp-ones(1,tp)*(yp')/tp;
        for s=1:kx
            xp=X(i,Tx(i,1,s):Tx(i,2,s),s);
            tp=Tx(i,2,s)-Tx(i,1,s)+1;
            Xd(i,Tx(i,1,s):Tx(i,2,s),s)=xp-ones(1,tp)*(xp')/tp;
        end %s
        for sz=1:kz
            for s=1:ktau
                xp=Z(i,Tz(i,1,sz,s):Tz(i,2,sz,s),sz,s);
                tp=Tz(i,2,sz,s)-Tz(i,1,sz,s)+1;
                Zd(i,Tz(i,1,sz,s):Tz(i,2,sz,s),sz,s)=xp-ones(1,tp)*(xp')/tp;
            end
        end %sz
    end
end

% create differences
dX=NaN(size(X));
dY=NaN(size(Y));
dYd=NaN(size(Y));
for s=1:kx
    dX(:,2:Tmax,s)=X(:,2:Tmax,s)-X(:,1:Tmax-1,s);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dY(:,2:Tmax)=Y(:,2:Tmax)-Y(:,1:Tmax-1);
%dY=Y;
dZ=NaN(N,Tmax,kz,ktau);
dZ(:,2:Tmax,:,:)=Z(:,2:Tmax,:,:)-Z(:,1:Tmax-1,:,:);

% compute star variables
dys=zeros(1,Tmax);
xs=zeros(k,Tmax);
zs=zeros(ktau,Tmax,kz);
dYn=dY; dYn(isnan(dY))=0;
Xdn=Xd; Xdn(isnan(Xd))=0;
Zdn=Zd; Zdn(isnan(Zd))=0;

for s=1:Tmax
   w=not(isnan(dY(:,s)));
   w=w/sum(w);
   dys(1,s)=w'*dYn(:,s);
end %s

for s=1:Tmax
   for h=1:k
      w=not(isnan(Xd(:,s,h)));
      w=w/sum(w);
      xs(h,s)=w'*Xdn(:,s,h);
   end %k
   for sz=1:kz
    for h=1:ktau
      w=not(isnan(Zd(:,s,sz,h)));
      w=w/sum(w);
      zs(h,s,sz)=w'*Zdn(:,s,sz,h);
    end %k
   end %sz
end %s



% largest lag
ml=max([lx, ldy, lcsx, lcsdy, lcsz, lz]);


i=countryi;

timin=max([max([Ty(i,1)+ldy+1,lcsdy+1]),max([max(Tx(i,1,:))+lx, lcsx+1]), max([ max(max(Tz(i,1,:,:)))+lz, lcsz+1])   ]);
timax=min([Ty(i,2),min(Tx(i,2,:)),  min(min(Tz(i,2,:,:))) ]);

    ti=timax-timin+1;
    tim=ti-1;
    % create contemporaneous star variables
    %ysc=ys(1,timin:timin+tim)';
    dysc=dys(1,timin:timin+tim)';
    for h=1:kx
        xsc(:,h)=xs(h,timin:timin+tim)';
    end
    for sz=1:kz
        for h=1:ktau
            zi(:,sz,h)=squeeze(Z(i,timin:timin+tim,sz,h))'; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        end
    end %sz    
    if lcsdy<0
        dysc=zeros(ti,0);
    end
    if lcsx<0
        xsc=zeros(ti,0);
    end
    
    
    % create lags for star variables
    dysl=zeros(ti, lcsdy);
    for s=1:lcsdy
        dysl(:,s)=dys(1,timin-s:timin-s+tim)';
    end
    xsl=zeros(ti, lcsx*kx);
    for s=1:lcsx
        for h=1:kx
            xsl(:,(s-1)*kx+h)=xs(h,timin-s:timin-s+tim)';
        end
    end
    zslii=zeros(ti, lcsz+1); 
    zsa=cell(kz,ktau);
   for sz=1:kz 
    for h=1:ktau
        for s=0:lcsz
            zslii(:,s+1)=zs(h,timin-s:timin-s+tim,sz)';
        end
        if lcsz>-1
            zsa{sz,h}=zslii;
        else
            zsa{sz,h}=zeros(ti,0);
        end
    end
   end %sz
    % create dependent variable
    dyi=dY(i,timin:timax)';
    
    % create regressors
        % first contemporaneous levels
    Xi=zeros(ti,kx);
    dyil=zeros(ti,ldy);
    for h=1:kx
        Xi(:,h)=X(i,timin:timax,h)';
    end
    if lx<0
       Xi=zeros(ti,0);
    end
    for s=1:ldy
        dyil(:,s)=dY(i,timin-s:timin-s+tim)';
    end
        % contemporaneous differences and lags
        xlags=zeros(ti, (lx)*kx);
        for s=1:lx
            for h=1:kx
               %xlags(:,(s-1)*kx+h)=X(i,timin-s:timin-s+tim, h)';
               % in differences instead
               xlags(:,(s-1)*kx+h)=X(i,timin-s+1:timin-s+tim+1, h)'-X(i,timin-s:timin-s+tim, h)';
            end% kx
        end %lx
     %create zlags   
    zlagslii=zeros(ti, lz); 
    zlagsa=cell(kz,ktau);
   for sz=1:kz 
    for h=1:ktau
        if dflz==1 % differences
            for s=1:lz
                zlagslii(:,s)=dZ(i,timin-s+1:timin-s+tim+1,sz, h)';  % Z...(N,Tmax,kz,ktau);
            end
        else % levels
            for s=1:lz
                zlagslii(:,s)=Z(i,timin-s:timin-s+tim,sz, h)';
            end   
        end
        zlagsa{sz,h}=zlagslii;
    end  
   end %sz        
        
    % create constant
    tau=ones(ti,1);
 
    XX=[ tau dyil Xi xlags dysc xsc dysl xsl ]; %  xs(1:T-1) xs(2:T)
    if size(Xi,2)>0 % Xi is present
        pos_xi=[size(tau,2)+size(dyil,2)+1,size(tau,2)+size(dyil,2)+size(Xi,2)];
    else
        pos_xi=[0,0];
    end
    if size(dyil,2)>0 % dyil is present
        pos_ylags=[size(tau,2)+1,size(tau,2)+size(dyil,2)];
    else
        pos_ylags=[0,0];
    end  
    
    Stars=[ dysc xsc dysl xsl ];
    if size(Stars,2)==0
        pos_stars=[0,0];
    else
        ss=size([tau dyil Xi xlags],2)+1;
        se=size(XX,2);
        pos_stars=[ss,se];
    end
    
xi=XX;
%debti=Xi(:,2);

















